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Abstract 

A simple model for the solidification of globular grains in metallic alloys is presented. 
Based on the Voronoi diagram of the nuclei centers, it accounts for the curvature of 
the grains near triple junctions. The predictions of this model are close to those of 
more refined approaches such as the phase field method, but with a computation cost 
decreased by several orders of magnitude. Therefore, this model is ideally suited for 
granular simulations linking the behavior of individual grains to macroscopic properties 
of the material. 

1 Introduction 

Microstructures that form during solidification of metallic alloys play a key role for the final 
properties of as-cast materials and for subsequent heat treatments. They also condition 
the formation of defects such as porosity and hot cracking [1]. As solidification of one or 
more solid phases is a moving free-boundary problem, several techniques have been developed 
especially over the past decade to overcome the difficulty of tracking the solid-liquid interface. 

One of the most powerful and widely used methods is the phase field technique [2] . In this 
method, the sharp interface is replaced by a continuous field varying from in the solid to 1 in 
the liquid with a diffuse interface over a finite thickness. As the phase field method is based 
on thermodynamics considerations, it provides a unified framework for many phenomena 
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(formation of dendrites, eutectics, peritectics, etc.). Yet, this method is very computation 
intensive as the mesh size has to be small with respect to the diffuse interface thickness, 
and the diffuse interface thickness has to be small with respect to the typical radius of 
curvature of the microstructure (e.g., dendrite tip radius, eutectic spacing). Moreover, since 
the problem is generally solved using an explicit scheme, a Fourier criterion must be satisfied 
to ensure numerical stability. This is why nowadays, phase field simulations are limited to the 
simulation of a few grains, very often in two dimensions. Other techniques such as the level 
set [3] or pseudo-front tracking [4] methods have been developed, with different advantages 
and drawbacks, but they roughly require the same amount of computations. 

If fundamental aspects of microstructure formation can be tackled with such sophisticated 
methods, it must also be recognized that there is an increasing need for simulations that can 
relate microstructure and macroscopic properties. This is particularly important for the last 
stage of solidification, during which the most important defects (porosity, hot tears) form. 
In particular, the gradual topological transition of the microstructure from continuous liquid 
films to a continuous and fully coherent solid is essential for hot tearing. Because of the 
random nature of nucleation, this transition is associated with the formation of increasingly 
larger clusters of grains as described by the percolation theory [5] . Therefore, the mechanical 
properties and feeding ability of the mushy zone during this transition cannot be deduced 
from the morphology of just a few grains [6]. 

In order to study the gradual formation of a coherent solid phase in globular alloys, 
an original approach has been proposed by Mathier et al. [7]. In this model, grains are 
approximated by polyhedrons based on the Voronoi diagram of a random set of nuclei. In 
order to compute the solidification of large and non-isothermal mushy zones, Vernede et 
al. further simphfied the assumptions of this model [8]. Yet, these two approaches lead 
to polyhedral grains and cannot account for the formation of isolated liquid pockets at the 
triple junctions of the grains. 

In the present contribution, a solidification model based on the approach of Mathier et al. 
but accounting for the Gibbs- Thomson effect near grain corners is described. The predictions 
of this model are shown to be close to those obtained with more refined methods such as the 
phase field or pseudo-front tracking methods, but with a computation time several orders 
of magnitude lower. After a brief recall of the granular m,odel of solidification based on 
the Voronoi tessellation, the account of the Gibbs-Thomson effect and its implementation 
in the model are described. In the last section, simulation results are compared with those 
obtained with a pseudo-front tracking method, thus providing a basis of discussion for setting 
up bonds to this simple approach. 

2 A solidification model based on Voronoi diagrams 

This 2 dimensional solidification model has been derived first by Mathier et al. [7] and 
further developed by Vernede et al. [8] . For a better understanding the main lines of this 
model are recalled here. 

The model assumes simultaneous nucleation of grains in a plane with a given density of 
random sites. Further assuming that the temperature difference across the average grain 
size is small with respect to the undercooling (i.e., small thermal gradient), the final grain 
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Figure 1: Various enlargements of the granular model: Voronoi tessellation associated with 
the nuclei centers (a); shape of the grains during solidification (b); solute balance within one 
triangle (c); smoothing procedure of the solid-liquid interface near the grain corners (d). 



structure is close to the Voronoi tessellation of the set of nuclei (Fig. 1, (a)) [9]. In the 
present work, the Voronoi tessellation was computed using the free access software qhull 
[10]. 

In order to further simplify the solidification model, the solute flux between elementary 
triangles is neglected in a first step. Thus, the smooth interface of each grain during growth 
can be approximated by a linear segment in each triangle connecting the nucleation center 
with a Voronoi segment. By construction, these segments are perpendicular to the vectors 
connecting the nucleation centers and the two triangles issued from the same Voronoi segment 
are symmetrical, (Fig. 1, (b)). 

Solidification is reduced therefore to a one-dimensional problem in each triangle, with the 
assumption of complete mixing of solute in the liquid phase and back-diffusion in the solid. 
Moreover, the temperature of the system is imposed, either uniform or given by a fixed 
thermal gradient and decreasing with a given cooling rate. Therefore, the solute balance 
integrated over the liquid phase of an elementary triangle gives: 

x*Ds^ix*) + v*x*{k - 1)q + ^{L' -x*')^ = (1) 

where x* and v* are the position and speed of the interface, respectively, Cg and q the solute 
concentration in the solid and liquid phase, c* the solute concentration in the solid at the 
interface, k the partition coefficient, t the time, Dg the diffusion coefficient in the solid, L the 
height of the elementary triangle perpendicular to the Voronoi segment and x the coordinate 
along this direction (Fig. 1, (c)). The last term of Eq. 1 accounts for the evolution of solute 
concentration, which is imposed by an external cooling rate T and the phase diagram. Thus 

x*Dg^{x*) + v*x*{k - 1)q +1(1^- X*")- = (2) 
ox ' ' 2 m 

where m is the slope of the liquidus. 

The first term in Eq. 2 associated with back-diffusion in the solid is computed by solving 

the cylindrical diffusion equation in the solid phase, which can be viewed as the solute balance 

between slices of the elementary triangle for constant x: 

dcs ^ ( d'^Cs 1 dcs \ , . 

dt ^ \ dx"^ X dx J 
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In order to easily account for solidification, we use a Landau transformation of the solid 
domain [0,x*(t)] into the reference ID domain [0,1], as introduced by VoUer and Sundarraj 
[11]: 



X 



Cs{x,t) ^ Cs{r),t) 77=— (4) 

X* 



Thus Eq. 3 becomes: 



dt J x*"^ drj 




(5) 



where the term accounts for the advection of the mesh with the sohdification front. This 
equation is solved using a finite difference scheme, with a nil fiux condition for 77 = and an 
imposed concentration c* — kci for = 1. Knowing the flux associated with back-diffusion 
and the concentration evolution in the liquid, the second term of Eq. 2 allows to deduce the 
velocity of the interface, v*, and thus to find the new position x*{t + dt). 

Note that solidification does not depend on the opening of the elementary triangles but 
only on its height L, i.e, on the half-distance between two nucleation centers. Near the end of 
solidification, the excess free energy needed to form a grain boundary out of two misoriented 
grains, i.e., coalescence undercooling [12], can be accounted for. This procedure is presented 
in Refs [7, 8] , but as this feature is not useful for the present contribution, it will not be 
detailed here. 



3 A model for grain corners 

Although the previous model predicts fairly well the evolutions of the grains, of the solid- 
liquid interface and of the solid fraction (see section 4), it leads to polyhedral grains (Fig. 1, 
(d)). This situation is clearly unrealistic as sharp corners are normally remelted by surface 
tension in non-faceted crystals. Furthermore, no liquid pocket can form at the triple junctions 
of the grains. As a consequence, this model overestimates the volume fraction of solid at 
which contact between grains occurs, as compared with experiments [8, 13]. 

In order to remove these hmitations, the Gibbs-Thomson effect at grain corners has to 
be considered. Assuming again that temperature is homogeneous at the scale of a grain, the 
liquid concentration at an interface with a local radius of curvature R is given by: 

where is the liquid solute concentration for a flat interface and F^^ is the Gibbs-Thomson 
coefficient. (Note that, in general, the liquidus slope m is negative). If solidiflcation would 
be arrested, solute will flow from low curvature areas to high curvature for coars- 

ening, thus remelting the highly curved zones. Yet, globular grains arc not spherical during 
solidiflcation and another phenomenon should balance the Gibbs-Thomson effect. 

Therefore, the idealized situation represented in Fig. 2 is considered for smoothing the 
shape of polyhedral grains based on the Gibbs-Thomson effect. The grain corner is modelled 
by a curved interface with a constant radius of curvature R, whereas elsewhere the interface 
is supposed flat and parallel to the flnal grain boundary. Moreover, the liquid is divided 
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into two zones, delimited by the dotted line AA' in Fig. 2. This hne, passes through the 
point separating the flat and circular portions of the solid-liquid interface (point labelled A 
in Fig. 2). Close to the corner, the hquid has a concentration equal to cf , whereas in the 
second zone, surrounding the planar interface, a homogeneous concentration equal to is 
considered. 

The gradient of solute around the corner can be estimated by (cf — c^)/(i? sin «) where 
i?sina is the distance between the tips of the corner and limit of the two zones. Moreover, 
we consider that the flux of solute is effective over a length of the order of R sin a. This leads 
to an estimation of the overall flux flowing from the zone surrounding the flat interface to 
the zone surrounding the corner: 

$~i?,(cr-cf) = -^>0 (7) 

vixtx 

We further look for a radius of curvature that is stable with time , i.e., i? = 0. This implies 
that the limit between the curved and planar interface stays along the dotted line AA' in 
Fig 2. Once this constrain is fixed, it can be shown that the rejection of solute associated 
with a fiat interface moving at a velocity is equal to that rejected by the projection of 
the curved interface moving at a velocity v^. Indeed, the projection of the curved interface 
is given by it! sin a, whereas the velocity is given by — v"^ / cos a. So one has: 

v°°R tan q; = i? sin av^ — Rsma (8) 

cos a 

where i?tana is the length of the extended fiat interface (dotted line AB in Fig. 2). How- 
ever, the grey surface is an additional volume of liquid, compared to the fiat interface, and 
represents the advantage of the corner for solute diffusion. In other words, the small incom- 
ing flux of solute contributes to increasing the concentration of the liquid in the area, 5', of 
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the grey zone. A more rigorous development, based on solute balances, is given in Appendix. 
But in summary, one has: 

$ = ^5 = dc^^ ^J_S^ l^(tana - a) (9) 
dt dt m m 2 ^ ' ^ ' 

Combining Eq. 7 and 9, one finally gets: 

R^^Ac ^ (10) 

tana — a —T 

where Ac is a dimensionless constant, arising from the simple description of the solute 
distribution (see Eq. 7 ). Nonetheless, simply setting its value to 1 produce satisfying results 
and this value will be used hereafter. 

It is interesting to note at this stage that Eq. 10 is close to a coarsening law: the radius 
of curvature of a grain corner is proportional to the third power of a driving force given 
by , where ATg is the solidification interval of the alloy, and to the third power of the 
solidification time, tf. The geometrical factor in front of this term is such that the radius of 
curvature becomes infinite when a is equal to (the grain corner is fiat) and nil when a is 
equal to | (the grain has disappeared). 

Using this relationship for the radius of curvature of the grain corners, the shape of the 
solid-liquid interface can be computed as follows. 

• The position of the flat interface or solid fraction is computed for each elementary 
triangle using the back-diffusion model described above 

• The radius of curvature at each grain corner is computed using Eq. 10. Please note 
that, for a fixed cooling rate, the radius is constant and can be computed only once 
before the time stepping procedure. 

• The interface in each elementary triangle is approximated by a flat portion and two 
rounded corners. If the length of the fiat interface becomes negative, the interface is 
approximated by an arc of a circle. 

• Rounded interfaces increase the overall volume of liquid by creating liquid pockets at 
grain corners. The position of the flat interface is then moved slightly forward in order 
to conserve the solid fraction computed with the fiat interface method (see Fig. 1 (d)). 



4 Validity and limits of the model 

The predictions of the present model have been compared with those obtained with a pseudo 
front tracking method (PFT) [4] . In this technique, the fraction of solid within each cell of 
an hexagonal network is computed based on an explicit solute diffusion calculation. A layer 
of cells always separate the solid and liquid phases and the solute fiux balance for such cells is 
converted into a solid fraction evolution. The position of the interface within these interfacial 
cells is computed using a piece-wise linear interface calculation (FLIC) algorithm [14]. Once 
the interface position within each interfacial cell is known, its curvature and the associated 
Gibss-Thomson effect arc calculated using a distance- field method similar to level-set. This 
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Figure 3: Solid fraction as a function of time computed with the Voronoi and the PFT 
methods 

method leads to predictions close to the phase field method as shown in more details in Ref. 



The same simulations have been carried out with the PFT and the Voronoi methods. 
Six nuclei have been randomly placed in a 4.5 ■ 10~^ mm^ domain with periodic boundary 
conditions. The average grain size is therefore around 90 ^m. An Al-lwt%Cu has been 
considered with a linearised phase diagram (m = —6.67 K/wt% and k = 0.14). The other 
parameters used in these calculations are: Ds = 1.5 • 10"^'^ m^/s, Di = 10~^ ■ m^/s, Tge = 
5 • 10-^ Km, T = -1 K/s. 

Figure 3 shows the evolution of the solid fraction, gs, computed with the two methods. 
Please keep in mind that, in the Voronoi method, the solid fraction calculated with rounded 
grains is equal to that obtained with polygonal grains. At the very beginning of solidification, 
grain growth predicted with the PFT method is slightly faster as the Voronoi method assumes 
complete mixing of solute in the liquid. As soon as a state of complete mixing is reached 
with the PFT method, the predictions of the two models are in very good agreement, even 
at high solid fraction. This shows that, despite the fairly strong assumptions of the (flat 
interface) Voronoi model, solute back-diffusion is well approximated. 

In order to compare in more details the shape of the interfaces predicted by the two 
methods, the normalized specific solid hquid interface, 5"°, is represented in Fig 4 as a 
function of Qs- This dimensionless number corresponds to the total length of the solid-liquid 
interface of the grains divided by the number of grains and by the average final grain size dc- 
This important parameter strongly influences the permeability of the grain assembly via the 
Carman-Kozeny relationship [15] as verified experimentally in aluminium alloys by Nielsen 
et al. [16]. A numerical study of S° also shows that this parameter is an indicator of the 
morphological transitions of the mushy zone [6] . As solid grains grow, the length of the solid 
liquid interface increases as in 2D until impingement/contact of the solid grains makes it 
go to zero. Thus, S° is maximum when the new contacts between the grains counterbalance 
the natural increase of the interface length. 

As can be seen on the left of Fig. 4, the overall shapes of S°{gs) are in excellent agreement. 



[4]. 



7 
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Figure 4: Comparison of the solid liquid interface shape predicted by the Voronoi and the 
PFT methods at various solid fractions. 



but the Voronoi method slightly overestimates this parameter, thus revealing that grains are 
still shghtly less rounded than those calculated with the PFT method. It should be specified 
that, under these conditions, the grains are clearly globular without significant destabilization 
of the interface. Formation of dendritic or globular dendritic grains would definitely increase 
the length of the solid-liquid interface, thus making S^ ppj^ > 5'° ^^,.. 

The grain shapes predicted by the two models at various solid fractions are also repre- 
sented in Fig. 4. For visualization purpose, the grains computed by the Voronoi method are 
represented in grey, whereas the interfaces predicted by the PFT method are represented 
with black lines. Again, a fairly good agreement between the two simulations, especially 
near the grain corners, can be seen. Please note that interfaces predicted with the Voronoi 
model are not necessarily continuous as solidification is computed separately for each elemen- 
tary triangle. Yet, these slight discontinuities do not affect much the topology of the liquid 
channel network, nor the estimation of the channels permeability and mushy zone topology. 

The maximum of 5*° is predicted to occur at Qs = 0.86 by the PFT model, while it is 
delayed to Qs = 0.89 with the Voronoi model. This difference can be easily understood by 
looking at the shape of the grains at Qs = 0.88. With the PFT model, the liquid film in 
between two close neighboring grains can break down into small droplets by coalescence, 
whereas with the Voronoi model it remains a film until final impingement. Such a liquid 
film instability occurs when the grain boundary energy is lower than twice the solid-liquid 
interfacial energy, i.e., attractive boundaries, and the film thickness is on the order of the 
thickness of the diffuse sohd-liquid interfaces. [12] No grain boundary energy has been set 
up in the PFT and Voronoi calculations, but the instability occurs too early with the PFT 
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Figure 5: Comparison of grain shape at (7^=0.68 for various sets of solidification parameters. 
Starting from the reference simulation (right figure), one parameter is changed at a time: 
Dg^ (a), Vge. (b) and k (c) are multiplied (top) or divided (bottom) by 2 with respect to the 
reference. 



method as it is conditioned by the mesh size rather than the diffuse interface thickness. Note 
that the quantitative simulation of the coalescence phenomenon remains a challenging topic 
despite advances done in the thermodynamic framework [12] and in the modelling of triple 
phase boundaries [17]. Please note also that grain boundary energies, and thus coalescence 
undercooling of repulsive interfaces, can be easily introduced in the present simulations. [8] 

Despite the numerical difficulties associated with simulation of the last stage solidification, 
the radii of CTirvature at grain corners calculated with Eq. 10 are close to those arising from 
the complex PFT calculation, even for relatively narrow channels. After breakdown of the 
liquid films, the sizes of the liquid pockets at triple junction points are also in fairly good 
agreement. Note that the isolated liquid pocket are represented nevertheless with a negative 
curvature in the Voronoi simulation. It is clear that the rounding procedure does not account 
for the coalescence phenomenon which transforms globally convex globular grains into convex 
liquid pockets. This simplified approach nevertheless allows to introduce realistic volumes 
of liquid at triple junctions, thus decreasing the volume fraction of solid at which the flat 
portions of the grain interfaces impinge. 

Figure 5 shows the effect of various parameters on the shape of the solid liquid interfaces. 
It is convenient to introduce a dimensionless number, C: 

^-^y (11) 

This number is the ratio of a "typical" curvature radius at grains corner (i.e., when 
(tana — a)^^ = 0.5) and the average grain size. Figures 5 a) and b) show that a larger 
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C number corresponds indeed to a larger radius of curvature at grains corners, regardless 
whether this is achieved by increasing D( or Vf^(. Moreover, the agreement between the PFT 
and Voronoi predictions remains when the physical parameters are changed. Yet, Fig. 5 also 
points out the limits of the present model, in particular at small C number. Indeed, the 
PFT predicts a destabilization of the globular grains into dendritic ones when C decreases. 
A criterion for the transition from globular to dendritic equiaxed grains has been recently 
proposed by Diepers et al.[18] These authors found that it occurs when: 



where Cp and L are the volumetric specific heat and latent of fusion, respectively, and 
Agd a dimensionless factor function of the anisotropy of the solid-liquid interfacial energy. 
Although this criterion has been derived for a sphere and a cooling rate imposed only at 
the boundaries of the system, it is interesting to note that the same power-law of the ratio 
^^^Y^ is retrieved. The additional term,(L/(A;AToCp)), comes from the solute undercooling 
ahead the solidification front and from an overall thermal balance. The destabilisation of the 
interface occurs at the beginning of solidification when this undercooling is crucial, whereas 
the selection of the grains corner curvature occurs latter in solidification when the solutal 
profile is fiat. A change in the partition factor /c, which does not affect the radius of curvature 
of the grains calculated with Eq. 10, nevertheless slightly infiuences the destabilization of the 
grains as predicted by Eq. 12 from the factor kHsTg and as observed with the PFT method 



For a high C number, another limitation of the model is encountered as illustrated by 
the solidification of the smallest grain in Figs. 5 a) b). The Voronoi model does not predict 
well the solidification of this grain at high C number. This is because the overall curvature 
of the grain is not accounted for in the Voronoi method. Yet, the curvature undercooling 
during growth of globulitic grains remains very small, typically 0.01 K for a spherical grain 
of 20 jim radius. That of the last fiquid droplets is slightly higher considering that the 
droplets are smaller and have a negative curvature, opposite to that considered here for the 
last liquid located at triple junctions. But in any case, curvature undercooling, despite its 
importance for coarsening and for the shape of the dendrite tip, is small in dendritic and 
globulitic solidification under normal conditions. 

In order to use a Voronoi construction, e.g., linear grain boundaries, to predict the final 
grain structures, two assumptions have been made: first, the temperature is uniform at the 
scale of the grains, and second, the nuclei all start at the same time. The first hypothesis is 
usually verified for globulitic structures, e.g., grains of about 100 /xm growing in a thermal 
gradient lower than 1 K/cm. The temperature difference across a reference specimen con- 
taining 100 grains is just 1 K in this case. The second hypothesis clearly pertains to alloys 
which are inoculated. It has been shown in the case of Al-alloys inoculated with TiB2 parti- 
cles that nucleation is athermal [19]. In athermal nucleation, the activation of a nucleant is 
a function of the undercooling only and not of time. In this case, the range of particle sizes 
which have been shown to be activated during solidification is typically 2-5 /xm in diameter 
[19]. Therefore, the maximum nucleation undercooling is about 0.2 K for a Gibbs- Thomson 
coefficient of 10~^ Km. This value is small compared to the growth undercoofing and thus 




(12) 



(Fig. 5 c). 
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an instantaneous nucleation assumption appears reasonable. 

The most striking feature of the present Voronoi model is its low computation cost. The 
PFT simulations presented in Fig. 4 take around 12 Hours on a 2.8 GHz Pentium 4 personal 
computer, whereas the Voronoi calculations requires less than 2 s! This represents a gain of 
more than 4 orders of magnitude. This difference is even more striking with large mushy 
zones. The computation of a whole mushy zone solidification that contains 14000 grains (see 
Fig. 9 of Ref. [6]) requires less than 10 seconds on the same computer. 



5 Conclusion 

A model for the solidification of globular grains based on Voronoi diagrams and round corners 
has been derived. This model is shown to estimate well the shape of the globular grains during 
solidification and to take into account the effect of the various physical parameters, providing 
the interface does not become unstable. Its very low computation cost makes it an ideal 

base for granular simulation of mushy zone at a mesoscale, that can compute macroscopic 
properties such as mechanical properties or feeding properties based on the behaviour of 
individual grains [6]. Therefore, this model can be easily extended to 3 dimensions, providing 
new opportunities for realistic granular simulation. 
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7 Appendix 

Considering a section of the flat interface far from the corner, the speed of the interface can 
be estimated from the following flux balance. 

h^ = cT{l-k)v--jra (13) 

The left hand term represents the variation of solute in the liquid, where h is the thickness 
of the liquid film. The first right hand term correspond to the solute rejected in the liquid 
due to the advance of the interface and is the flux pumped in the solid by back-diffusion. 
Similarly, a solute balance on the liquid part surrounding the corner (see main section) can 
be derived. If Q denotes this domain, one has. 



/ 

JQ 



^ds = cf{l-k)Rsmav^ -aRj^a + ^ (14) 

Please note that this balance accounts for the flux $ exchanged between the flat and curved 
portions of the interface. 
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The variation of solute concentration at the interface is imposed by the coohng rate and 
dtcf^ = dfcf — T jm. As a consequence, whatever is the precise repartition of solute around 
the grain corner, the variation of solute around the grain corner can be estimated by: 

/ ^-^ds = Sn- (15) 
JQ at m 

where is the area of the domain surrounding the round corner. As v'^ — v'^ j cos a (see 
main section), Eqs 13, 14 and 15 give: 

Sn- = ^Rtana{h-+jr,) - aRj^, + ^ (16) 

m cf m 

Considering that cf/cf ~ 1 and neglecting the differences of back diffusion along the flat 
and curved parts of the interface, one gets: 

T T 

(S + Rh tana)— = Rh tan a h$ (17) 

m m 

Where the surface Sq have been separated into the surface dehmited by the extension of 
the flat interface {Rhtajaa) and an extra surface S represented in grey in Fig. 2. As stated 
in the main part, one retrieve the fact that the solute rejected by the flat interface moving 
at velocity is equivalent to that rejected by the curved interface moving at velocity 
yT _ I (,QgQ, Removing this term on the left and right hand sides finally gives: 

5- = $ (18) 
m 

This represents the solute balance between the solute flux induced by the Gibbs- Thomson 
effect and the geometrical advantage of a corner for diffusion. 
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